####
#    This syntax replicates the results of the paper
#    Gr�mping, Max. 2019. "More Bang for the Buck: Media Freedom and Organizational Strategies in the Agenda-Setting of Human Rights Groups",
#    in: Political Communication, doi: 10.1080/10584609.2018.1551256.
# 
#    The dataset includes five variables that are additive indices, for which missing data was imputed: 
#    ("DEMI.EFFORT.i", "DEMI.FOCUS.i", "DEMI.PROFESS.i", "DEMI.RESOURCES.i", "DEMI.EXPERIENCE.i")
#    The multiple imputation procedure and creation of indices is detailed in the separate syntax "Imputation of IVs.R",
#    based on raw survey responses in "Indata-Imputation of IVs.RData"
####

rm(list=ls())

require(pacman)
pacman::p_load(dplyr)
pacman::p_load(Hmisc)
pacman::p_load(ggplot2)
pacman::p_load(glmmTMB)
pacman::p_load(lme4)
pacman::p_load(optimx)
pacman::p_load(moments)
pacman::p_load(nsRFA)
pacman::p_load(texreg)
pacman::p_load(bbmle)
pacman::p_load(reshape2)
pacman::p_load(broom)
pacman::p_load(stringr)
pacman::p_load(lemon)
pacman::p_load(merTools)
pacman::p_load(polycor)
pacman::p_load(psych)
pacman::p_load(stargazer)
pacman::p_load(sjPlot)
pacman::p_load(dataMaid)


theme_set(theme_bw())
options(scipen=999)

setwd("C:\\ENTER PATH TO WORKING DIRECTORY\\")

load("Replication data-UPCP-2018-0112.RData")

# Create functions --------------------------------------------------------


# Multiple plot function   (http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/)
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}



# Function that returns Root Mean Squared Error
rmse <- function(error)
{
  sqrt(mean(error^2))
}



##############  ggplot decimals display function
fmt_dcimals <- function(decimals=0){
  # return a function responpsible for formatting the 
  # axis labels with a given number of decimals 
  function(x) as.character(round(x,decimals))
}





# Variables ----------------------------------------------------------------

#   	ID	=	Unique ID of group
#   	ISO	=	ISO country code
#   	election	=	Election for which attention is measured (ISO_DDMMYYYY_Presidential/Legislative)
#   	att_dom	=	Domestic news attention, raw count
#   	att_dom_std	=	Domestic news attention, standardized
#   	att_dom_bin	=	Domestic news attention, dichotomized
#   	att_intl         	=	Int'l news attention, raw count
#   	att_intl_std	=	Int'l news attention, standardized
#   	att_intl_bin	=	Int'l news attention, dichotomized
#   	rsf_pfi	=	Press freedom (RSF)
#   	rating.inv	=	Issue urgency (extent of electoral malpractice) (PEI)
#   	DEMI.EFFORT.i	=	Media effort, imputed
#   	DEMI.FOCUS.i	=	Specialization, imputed
#   	DEMI.PROFESS.i   	=	Professionalization, imputed
#   	DEMI.EFFORT	=	Media effort, raw
#   	DEMI.FOCUS	=	Specialization, raw
#   	DEMI.PROFESS	=	Professionalization, raw
#   	DEMI.RESOURCES.i	=	Resources
#   	DEMI.EXPERIENCE.i	=	Longevity
#   	DEMI.1	=	Network membership
#   	att_dom_bin_lag  	=	Lagged domestic news attention, dichotomized
#   	att_intl_bin_lag	=	Lagged int'l news attention, dichotomized
#   	att_dom_lag	=	Lagged domestic news attention, raw count
#   	att_intl_lag	=	Lagged int'l news attention, raw count
#   	att_elect_dom	=	Domestic attention to elections (offset in M3)
#   	att_elect_intl	=	Int'l attention to elections (offset in M4)
#   	fhcategory	=	Freedom House status (FH)
#   	PEItype          	=	Election quality (PEI)
#   	unna_gdppc	=	GDP per capita (PPP) in adjusted 2009 US$ (UN)
#   	undp_hdi	=	Human Development Index (UNDP)
#   	wdi_internet	=	Internet penetration (World Bank)
#   	p_polity2	=	Revised Combined Polity Score
#   	unna_pop	=	Population (UN)
#   	A37_worked	=	Time worked for organization (1 = <1yr; 2 = 1-2yrs; 3 = 3-5yrs; 4 = 6-10Yrs: 5 = 11-20yrs; 6 = >20yrs)
#   	A38_strategic    	=	Involvement in strategic decision (0= not at all involved; 10 = very involved)
#   	A39_operations	=	Involvement in operational decisions (0= not at all involved; 10 = very involved)
#   	att_dom2	=	Alternative measure for domestic news attention, raw count
#   	att_dom2_std 	=	Alternative measure for domestic news attention, standardized


# Models M1: explaining the threshold of domestic attention to groups  --------
M1null <- glmer(att_dom_bin ~ 1 + (1 | ISO), data=df, family=binomial(link = "logit"),
                     control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
M1a <- glmer(att_dom_bin ~ 
                     scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) + 
                     (1 | ISO), data=df, family=binomial(link = "logit"),
                     control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
M1b <- glmer(att_dom_bin ~ 
                     scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) + 
                     scale(att_dom_bin_lag) + 
                     (1 | ISO), data=df, family=binomial(link = "logit"),
                     control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
M1c <- glmer(att_dom_bin ~ 
                     scale(rsf_pfi)*scale(rating.inv) + 
                     scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) + 
                     scale(att_dom_bin_lag) +
                     (1 | ISO), data=df, family=binomial(link = "logit"),
                     control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
M1d <- glmer(att_dom_bin ~ scale(DEMI.EFFORT.i)  + scale(DEMI.FOCUS.i) + scale(DEMI.PROFESS.i) + 
                     scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) + 
                     scale(att_dom_bin_lag) +
                     (1 | ISO), data=df, family=binomial(link = "logit"),
                     control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
M1e <- glmer(att_dom_bin ~ scale(rsf_pfi)*scale(rating.inv) + 
                     scale(DEMI.EFFORT.i)  + scale(DEMI.FOCUS.i) + scale(DEMI.PROFESS.i) +
                     scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) + 
                     scale(att_dom_bin_lag) +
                     (1 | ISO), data=df, family=binomial(link = "logit"),
                     control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
M1 <- glmer(att_dom_bin ~ scale(rsf_pfi)*scale(rating.inv) +
                    scale(rsf_pfi)*scale(DEMI.EFFORT.i)  + scale(rsf_pfi)*scale(DEMI.FOCUS.i) + scale(rsf_pfi)*scale(DEMI.PROFESS.i) + 
                    scale(DEMI.RESOURCES.i) +scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) +
                    scale(att_dom_bin_lag) +  
                    (1 | ISO), data=df, family=binomial(link = "logit"),
                    control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))

# Models M2: explaining the threshold of int'l attention to groups  --------
M2null <- glmer(att_intl_bin ~ 1 + (1 | ISO), data=df, family=binomial(link = "logit"),
                     control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
M2a <- glmer(att_intl_bin ~ 
                     scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) + 
                     (1 | ISO), data=df, family=binomial(link = "logit"),
                     control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
M2b <- glmer(att_intl_bin ~ 
                     scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) + 
                     scale(att_intl_bin_lag) + 
                     (1 | ISO), data=df, family=binomial(link = "logit"),
                     control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
M2c <- glmer(att_intl_bin ~ 
                     scale(rsf_pfi)*scale(rating.inv) + 
                     scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) + 
                     scale(att_intl_bin_lag) +
                     (1 | ISO), data=df, family=binomial(link = "logit"),
                     control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
M2d <- glmer(att_intl_bin ~ scale(DEMI.EFFORT.i)  + scale(DEMI.FOCUS.i) + scale(DEMI.PROFESS.i) + 
                     scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) + 
                     scale(att_intl_bin_lag) +
                     (1 | ISO), data=df, family=binomial(link = "logit"),
                     control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
M2e <- glmer(att_intl_bin ~ scale(rsf_pfi)*scale(rating.inv) + 
                     scale(DEMI.EFFORT.i)  + scale(DEMI.FOCUS.i) + scale(DEMI.PROFESS.i) +
                     scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) + 
                     scale(att_intl_bin_lag) +
                     (1 | ISO), data=df, family=binomial(link = "logit"),
                     control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
M2 <- glmer(att_intl_bin ~ scale(rsf_pfi)*scale(rating.inv) +
                    scale(rsf_pfi)*scale(DEMI.EFFORT.i)  + scale(rsf_pfi)*scale(DEMI.FOCUS.i) + scale(rsf_pfi)*scale(DEMI.PROFESS.i) + 
                    scale(DEMI.RESOURCES.i) +scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) +
                    scale(att_intl_bin_lag) +  
                    (1 | ISO), data=df, family=binomial(link = "logit"),
                    control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))

# Models M3: explaining the count of domestic attention to groups -------- 

dfc3 <- df[complete.cases(df[c("att_dom",
                                            "DEMI.RESOURCES.i",
                                            "DEMI.EFFORT.i",
                                            "DEMI.FOCUS.i",
                                            "DEMI.PROFESS.i",
                                            "DEMI.EXPERIENCE.i",
                                            "DEMI.1",
                                            "rsf_pfi",
                                            "rating.inv",
                                            "att_dom_lag",
                                            "p_polity2",
                                            "att_elect_dom",
                                            "unna_gdppc",
                                            "p_polity2",
                                            "ISO")]),]

form <- att_dom ~ scale(rsf_pfi)*scale(rating.inv) + 
  scale(rsf_pfi)*scale(DEMI.EFFORT.i)  + scale(rsf_pfi)*scale(DEMI.FOCUS.i) + scale(rsf_pfi)*scale(DEMI.PROFESS.i) + 
  scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) + 
  scale(att_dom_lag) + 
  offset(log(att_elect_dom)) +
  (1|ISO)
zi <- ~ scale(rsf_pfi)*scale(rating.inv) +
  scale(DEMI.EFFORT.i)  + scale(DEMI.FOCUS.i) +
  scale(att_dom_lag) + offset(log(att_elect_dom)) +
  (1|ISO)

M3p = glmmTMB(form, data=dfc3, family=poisson, control = glmmTMBControl(optCtrl=list(iter.max=1e9,eval.max=1e9)))
M3nb = glmmTMB(form, data=dfc3, family=nbinom2, control = glmmTMBControl(optCtrl=list(iter.max=1e9,eval.max=1e9)))
M3zinb = glmmTMB(form, zi=zi, dfc3, family=nbinom2, control = glmmTMBControl(optCtrl=list(iter.max=1e9,eval.max=1e9))) #does not converge
M3hp = glmmTMB(form, zi = zi, dfc3, family=truncated_poisson, control = glmmTMBControl(optCtrl=list(iter.max=1e9,eval.max=1e9)))
M3hnb = glmmTMB(form, zi = zi, dfc3, family=truncated_nbinom2, control = glmmTMBControl(optCtrl=list(iter.max=1e9,eval.max=1e9)))

bbmle::AICtab(M3p, M3nb, M3zinb, M3hp, M3hnb)

rm(form)
rm(zi)



# Models M4: explaining the count of int'l attention to groups  --------

dfc4 <- df[complete.cases(df[c("att_intl",
                               "DEMI.RESOURCES.i",
                               "DEMI.EFFORT.i",
                               "DEMI.FOCUS.i",
                               "DEMI.PROFESS.i",
                               "DEMI.EXPERIENCE.i",
                               "DEMI.1",
                               "rsf_pfi",
                               "rating.inv",
                               "att_intl_lag",
                               "att_elect_intl",
                               "unna_gdppc",
                               "p_polity2",
                               "ISO")]),]

form <- att_intl ~ scale(rsf_pfi)*scale(rating.inv) + 
  scale(rsf_pfi)*scale(DEMI.EFFORT.i)  + scale(rsf_pfi)*scale(DEMI.FOCUS.i) + scale(rsf_pfi)*scale(DEMI.PROFESS.i) + 
  scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) +
  scale(att_intl_lag) + 
  offset(log(att_elect_intl)) +
  (1|ISO)
zi <- ~ scale(rsf_pfi)*scale(rating.inv) +
  scale(DEMI.EFFORT.i)  + scale(DEMI.FOCUS.i) +
  scale(att_intl_lag) + offset(log(att_elect_intl)) +
  (1|ISO)

M4p = glmmTMB(form, data=dfc4, family=poisson, control = glmmTMBControl(optCtrl=list(iter.max=1e9,eval.max=1e9)))
M4nb = glmmTMB(form, data=dfc4, family=nbinom2, control = glmmTMBControl(optCtrl=list(iter.max=1e9,eval.max=1e9)))
M4hp = glmmTMB(form, zi = zi, dfc4, family=truncated_poisson, control = glmmTMBControl(optCtrl=list(iter.max=1e9,eval.max=1e9)))
M4hnb = glmmTMB(form, zi = zi, dfc4, family=truncated_nbinom2, control = glmmTMBControl(optCtrl=list(iter.max=1e9,eval.max=1e9)))
M4zinb = glmmTMB(form, zi=zi, dfc4, family=nbinom2, control = glmmTMBControl(optCtrl=list(iter.max=1e9,eval.max=1e9)))  #does not converge

bbmle::AICtab(M4p, M4nb, M4hp, M4hnb)

rm(form)
rm(zi)


# Figure 1: Cumulative distribution of media attention at group level --------

df_prop <- df[c("att_dom_std", "att_intl_std")]
df_prop <- df_prop[order(df_prop$att_intl_std),]
df_prop <- melt(df_prop)

ggplot(df_prop, aes(value, linetype = as.factor(variable))) +
  stat_ecdf(geom = "line", pad = FALSE) +
  xlab("Proportion of media attention") +
  ylab("Percent of groups") +
  scale_y_continuous(label=scales::percent_format(accuracy = 1)) +
  scale_x_continuous(label=scales::percent_format(accuracy = 1)) +
  scale_linetype_discrete(name="Arena",
                          breaks=c("att_dom_std", "att_intl_std"),
                          labels=c("Domestic news", "Int'l news")) +
  coord_cartesian(xlim = c(0,0.20)) 

ggsave(filename = "Fig 1.jpg", width = 5, height = 3, dpi = 600, units = "in")

rm(df_prop)





# Figure 2: Cumulative distribution of skew and kurtosis of media attention at country level --------

dft <- df[c("ISO", "att_dom_std")]
dft <- dft[!is.na(dft$att_dom_std),]
dft <- dft[order(dft$att_dom_std,decreasing=TRUE),] 
dft <- transform(dft,rank=ave(1:nrow(dft),ISO,
                            FUN=function(x) order(att_dom_std[x],decreasing=TRUE)))
dft <- dft[order(dft$ISO, dft$rank),] 
dft <- dcast(dft, ISO ~ rank, value.var = "att_dom_std")
test <-as.data.frame(apply(dft[,2:144], 1, skewness, na.rm =TRUE))
names(test)[1] <- "att_dom_std.skew"
dft <- cbind(dft, test)
test <-as.data.frame(apply(dft[,2:144], 1, kurtosis, na.rm =TRUE))
names(test)[1] <- "att_dom_std.kurtosis"
dft <- cbind(dft, test)
att_dom_std <- dft[c("ISO" , "att_dom_std.skew","att_dom_std.kurtosis")]
rm(dft)
rm(test)

dft <- df[c("ISO", "att_intl_std")]
dft <- dft[!is.na(dft$att_intl_std),]
dft <- dft[order(dft$att_intl_std,decreasing=TRUE),] 
dft <- transform(dft,rank=ave(1:nrow(dft),ISO,
                            FUN=function(x) order(att_intl_std[x],decreasing=TRUE)))
dft <- dft[order(dft$ISO, dft$rank),] 
dft <- dcast(dft, ISO ~ rank, value.var = "att_intl_std")
test <-as.data.frame(apply(dft[,2:144], 1, skewness, na.rm =TRUE))
names(test)[1] <- "att_intl_std.skew"
dft <- cbind(dft, test)
test <-as.data.frame(apply(dft[,2:144], 1, kurtosis, na.rm =TRUE))
names(test)[1] <- "att_intl_std.kurtosis"
dft <- cbind(dft, test)
att_intl_std <- dft[c("ISO" , "att_intl_std.skew","att_intl_std.kurtosis")]
rm(dft)
rm(test)

mom <- merge(att_dom_std, att_intl_std, by.x = "ISO", all=T)

sk <- mom[,c(2,4)]
sk <- melt(sk)

kur <- mom[,c(3,5)]
kur <- melt(kur)

g1 <- ggplot(sk, aes(value, linetype = as.factor(variable))) +
  stat_ecdf(geom = "line", pad = FALSE) +
  xlab("Skewness") +
  ylab("Percent of countries") +
  scale_y_continuous(label=scales::percent) +
  scale_linetype_discrete(name="Arena",
                          breaks=c("att_dom_std.skew", "att_intl_std.skew"),
                          labels=c("Domestic news", "Int'l news")) 
g2 <- ggplot(kur, aes(value, linetype = as.factor(variable))) +
  stat_ecdf(geom = "line", pad = FALSE) +
  xlab("Kurtosis") +
  ylab("Percent of countries") +
  scale_y_continuous(label=scales::percent) +
  scale_linetype_discrete(name="Arena",
                          breaks=c("att_dom_std.kurtosis", "att_intl_std.kurtosis"),
                          labels=c("Domestic news", "Int'l news")) 

ggsave(plot = multiplot(g1, g2), filename = "Fig 2.jpg", width = 5, height = 4, dpi = 600, units = "in")

rm(g1)
rm(g2)
rm(sk)
rm(kur)
rm(mom)
rm(att_dom_std)
rm(att_intl_std)






# Figure 3: Explaining media attention to human rights groups --------

# extract coefficients and name them
M1coef <- broom::tidy(M1)
M1coef <- M1coef[,c(1,2,3,5)]
colnames(M1coef) <- c("term", "estimate", "std.error", "p.value")
M1coef$model <- "M1"
M1coef$model2 <- "GLMM"
M1coef$term[M1coef$term=="(Intercept)"] <- "Intercept"
M1coef$term[M1coef$term=="scale(rating.inv)"] <- "H1a: Issue urgency"
M1coef$term[M1coef$term=="scale(rsf_pfi)"] <- "H1b: Press freedom"
M1coef$term[M1coef$term=="scale(rsf_pfi):scale(rating.inv)"] <- "H1c: Press freedom X Issue urgency"
M1coef$term[M1coef$term=="scale(DEMI.EFFORT.i)"] <- "H2a: Media effort"
M1coef$term[M1coef$term=="scale(DEMI.FOCUS.i)"] <- "H2b: Specialization"
M1coef$term[M1coef$term=="scale(DEMI.PROFESS.i)"] <- "H2c: Professionalization"
M1coef$term[M1coef$term=="scale(rsf_pfi):scale(DEMI.EFFORT.i)"] <- "H3a: Press freedom X Media effort"
M1coef$term[M1coef$term=="scale(rsf_pfi):scale(DEMI.FOCUS.i)"] <- "H3b: Press freedom X Specialization"
M1coef$term[M1coef$term=="scale(rsf_pfi):scale(DEMI.PROFESS.i)"] <- "H3c: Press freedom X Professionalization"
M1coef$term[M1coef$term=="scale(DEMI.RESOURCES.i)"] <- "Resources"
M1coef$term[M1coef$term=="scale(DEMI.EXPERIENCE.i)"] <- "Longevity"
M1coef$term[M1coef$term=="scale(DEMI.1)"] <- "Network membership"
M1coef$term[M1coef$term=="scale(att_dom_bin_lag)"] <- "Lagged media attention"

M2coef <- broom::tidy(M2)
M2coef <- M2coef[,c(1,2,3,5)]
colnames(M2coef) <- c("term", "estimate", "std.error", "p.value")
M2coef$model <- "M2"
M2coef$model2 <- "GLMM"
M2coef$term[M2coef$term=="(Intercept)"] <- "Intercept"
M2coef$term[M2coef$term=="scale(rating.inv)"] <- "H1a: Issue urgency"
M2coef$term[M2coef$term=="scale(rsf_pfi)"] <- "H1b: Press freedom"
M2coef$term[M2coef$term=="scale(rsf_pfi):scale(rating.inv)"] <- "H1c: Press freedom X Issue urgency"
M2coef$term[M2coef$term=="scale(DEMI.EFFORT.i)"] <- "H2a: Media effort"
M2coef$term[M2coef$term=="scale(DEMI.FOCUS.i)"] <- "H2b: Specialization"
M2coef$term[M2coef$term=="scale(DEMI.PROFESS.i)"] <- "H2c: Professionalization"
M2coef$term[M2coef$term=="scale(rsf_pfi):scale(DEMI.EFFORT.i)"] <- "H3a: Press freedom X Media effort"
M2coef$term[M2coef$term=="scale(rsf_pfi):scale(DEMI.FOCUS.i)"] <- "H3b: Press freedom X Specialization"
M2coef$term[M2coef$term=="scale(rsf_pfi):scale(DEMI.PROFESS.i)"] <- "H3c: Press freedom X Professionalization"
M2coef$term[M2coef$term=="scale(DEMI.RESOURCES.i)"] <- "Resources"
M2coef$term[M2coef$term=="scale(DEMI.EXPERIENCE.i)"] <- "Longevity"
M2coef$term[M2coef$term=="scale(DEMI.1)"] <- "Network membership"
M2coef$term[M2coef$term=="scale(att_intl_bin_lag)"] <- "Lagged media attention"

M3coef <- data.frame(term     = rep(names(coef(summary(M3nb))$cond[, "Estimate"])),
                     estimate = c(coef(summary(M3nb))$cond[, "Estimate"]),
                     std.error = c(coef(summary(M3nb))$cond[, "Std. Error"]),
                     p.value = c(coef(summary(M3nb))$cond[, "Pr(>|z|)"]),
                     model    = rep(c("M3"), each = 14)
)
M3coef$term <- as.character(M3coef$term)
M3coef$model2 <- "NB"
M3coef$term[M3coef$term=="(Intercept)"] <- "Intercept"
M3coef$term[M3coef$term=="scale(rating.inv)"] <- "H1a: Issue urgency"
M3coef$term[M3coef$term=="scale(rsf_pfi)"] <- "H1b: Press freedom"
M3coef$term[M3coef$term=="scale(rsf_pfi):scale(rating.inv)"] <- "H1c: Press freedom X Issue urgency"
M3coef$term[M3coef$term=="scale(DEMI.EFFORT.i)"] <- "H2a: Media effort"
M3coef$term[M3coef$term=="scale(DEMI.FOCUS.i)"] <- "H2b: Specialization"
M3coef$term[M3coef$term=="scale(DEMI.PROFESS.i)"] <- "H2c: Professionalization"
M3coef$term[M3coef$term=="scale(rsf_pfi):scale(DEMI.EFFORT.i)"] <- "H3a: Press freedom X Media effort"
M3coef$term[M3coef$term=="scale(rsf_pfi):scale(DEMI.FOCUS.i)"] <- "H3b: Press freedom X Specialization"
M3coef$term[M3coef$term=="scale(rsf_pfi):scale(DEMI.PROFESS.i)"] <- "H3c: Press freedom X Professionalization"
M3coef$term[M3coef$term=="scale(DEMI.RESOURCES.i)"] <- "Resources"
M3coef$term[M3coef$term=="scale(DEMI.EXPERIENCE.i)"] <- "Longevity"
M3coef$term[M3coef$term=="scale(DEMI.1)"] <- "Network membership"
M3coef$term[M3coef$term=="scale(att_dom_lag)"] <- "Lagged media attention"

M4coef <- data.frame(term     = rep(names(coef(summary(M4nb))$cond[, "Estimate"])),
                     estimate = c(coef(summary(M4nb))$cond[, "Estimate"]),
                     std.error = c(coef(summary(M4nb))$cond[, "Std. Error"]),
                     p.value = c(coef(summary(M4nb))$cond[, "Pr(>|z|)"]),
                     model    = rep(c("M4"), each = 14)
)
M4coef$term <- as.character(M4coef$term)
M4coef$model2 <- "NB"
M4coef$term[M4coef$term=="(Intercept)"] <- "Intercept"
M4coef$term[M4coef$term=="scale(rating.inv)"] <- "H1a: Issue urgency"
M4coef$term[M4coef$term=="scale(rsf_pfi)"] <- "H1b: Press freedom"
M4coef$term[M4coef$term=="scale(rsf_pfi):scale(rating.inv)"] <- "H1c: Press freedom X Issue urgency"
M4coef$term[M4coef$term=="scale(DEMI.EFFORT.i)"] <- "H2a: Media effort"
M4coef$term[M4coef$term=="scale(DEMI.FOCUS.i)"] <- "H2b: Specialization"
M4coef$term[M4coef$term=="scale(DEMI.PROFESS.i)"] <- "H2c: Professionalization"
M4coef$term[M4coef$term=="scale(rsf_pfi):scale(DEMI.EFFORT.i)"] <- "H3a: Press freedom X Media effort"
M4coef$term[M4coef$term=="scale(rsf_pfi):scale(DEMI.FOCUS.i)"] <- "H3b: Press freedom X Specialization"
M4coef$term[M4coef$term=="scale(rsf_pfi):scale(DEMI.PROFESS.i)"] <- "H3c: Press freedom X Professionalization"
M4coef$term[M4coef$term=="scale(DEMI.RESOURCES.i)"] <- "Resources"
M4coef$term[M4coef$term=="scale(DEMI.EXPERIENCE.i)"] <- "Longevity"
M4coef$term[M4coef$term=="scale(DEMI.1)"] <- "Network membership"
M4coef$term[M4coef$term=="scale(att_intl_lag)"] <- "Lagged media attention"

x <- c("H1a: Issue urgency",
       "H1b: Press freedom",
       "H1c: Press freedom X Issue urgency",
       "H2a: Media effort",
       "H2b: Specialization",
       "H2c: Professionalization",
       "H3a: Press freedom X Media effort",
       "H3b: Press freedom X Specialization",
       "H3c: Press freedom X Professionalization",
       "Resources",
       "Longevity",
       "Network membership",
       "Lagged media attention"
)

M1coef <- M1coef[match(x, M1coef$term),]
M2coef <- M2coef[match(x, M2coef$term),]
M3coef <- M3coef[match(x, M3coef$term),]
M4coef <- M4coef[match(x, M4coef$term),]

# Specify the width of confidence intervals
interval1 <- -qnorm((1-0.9)/2)  # 90% multiplier
interval2 <- -qnorm((1-0.95)/2)  # 95% multiplier

# Plot
finalcoefs <- rbind(M1coef, M2coef, M3coef, M4coef)
finalcoefs$term = factor(finalcoefs$term, levels=c("H1a: Issue urgency",
                                                   "H1b: Press freedom",
                                                   "H1c: Press freedom X Issue urgency",
                                                   "H2a: Media effort",
                                                   "H2b: Specialization",
                                                   "H2c: Professionalization",
                                                   "H3a: Press freedom X Media effort",
                                                   "H3b: Press freedom X Specialization",
                                                   "H3c: Press freedom X Professionalization",
                                                   "Resources",
                                                   "Longevity",
                                                   "Network membership",
                                                   "Lagged media attention"))
var_width = 30
finalcoefs <- mutate(finalcoefs, term.short = str_wrap(term, width = var_width))

ggplot(finalcoefs, aes(shape = model)) +
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) +
  geom_pointrange(aes(x = 0, y = estimate, ymin = estimate - std.error*interval2,
                      ymax = estimate + std.error*interval2),
                  position = position_dodge(width = 3),
                  fill = "WHITE") +
  scale_shape_manual(
    values = c(16,15,1,0),
    limits = c("M1", "M2", "M3", "M4"),
    labels= c("M1", "M2", "M3", "M4")
  ) +
  scale_x_discrete(breaks=NULL) +
  scale_y_continuous(breaks = c(-0.5, 0, 0.5, 1, 1.5)) +
  facet_rep_grid(term ~ ., switch = "x",  labeller = label_wrap_gen(width = 25, multi_line = TRUE))+
  theme_bw() + 
  theme(strip.text.y = element_text(angle = 0), strip.background = element_rect(colour="white", fill="white"))+
  theme(panel.border=element_blank(),
        axis.line = element_line(),
        axis.ticks = element_line(size = 1)) +
  coord_capped_cart(bottom='both', left='both') +
  coord_cartesian(ylim = c(-0.5, 1.5)) +
  coord_capped_flip()

ggsave(filename = "Fig 3.jpg", width = 5, height = 7, dpi = 600, units = "in")


rm(finalcoefs)
rm(M1coef)
rm(M2coef)
rm(M3coef)
rm(M4coef)
rm(interval1)
rm(interval2)
rm(x)
rm(var_width)


# Figure 4: The effect of institutional constraints and issue urgency on the predicted probability of achieving media attention --------

# predict
# see https://stackoverflow.com/questions/37015802/generating-marginal-prediction-confidence-intervals-from-a-glmer-object-using-pr
# see https://stackoverflow.com/questions/31948849/generating-predictive-simulations-from-a-multilevel-model-with-random-intercepts

# predict domestic attention
medEff = merTools::REquantile(M1, quantile = 0.5,groupFctr = "ISO", term = "(Intercept)")
newdata <- expand.grid(rating.inv=seq(4, 9, length.out=20),
                       rsf_pfi=c(0.4, 0.9),
                       DEMI.EFFORT.i=mean(df$DEMI.EFFORT.i, na.rm=T),
                       DEMI.FOCUS.i=mean(df$DEMI.FOCUS.i, na.rm=T),
                       DEMI.PROFESS.i=mean(df$DEMI.PROFESS.i, na.rm=T),
                       DEMI.RESOURCES.i=mean(df$DEMI.RESOURCES.i, na.rm=T),
                       DEMI.EXPERIENCE.i=mean(df$DEMI.EXPERIENCE.i, na.rm=T),
                       DEMI.1=mean(df$DEMI.1, na.rm=T),
                       att_dom_bin_lag=mean(df$att_dom_bin_lag, na.rm=T),
                       ISO=medEff)
out1 <- merTools::predictInterval(M1, newdata=newdata, level=0.90,
                                  n.sims=1000, stat="mean", type = "probability", include.resid.var=FALSE, fix.intercept.variance=T)
out1 <- cbind(out1, newdata$rating.inv, newdata$rsf_pfi)
names(out1) <- c("fit","upr","lwr","rating.inv","rsf_pfi" )
out1$arena <- "Domestic"

# predict int'l attention
medEff = merTools::REquantile(M2, quantile = 0.5,groupFctr = "ISO", term = "(Intercept)")
newdata <- expand.grid(rating.inv=seq(4, 9, length.out=20),
                       rsf_pfi=c(0.4, 0.9),
                       DEMI.EFFORT.i=mean(df$DEMI.EFFORT.i, na.rm=T),
                       DEMI.FOCUS.i=mean(df$DEMI.FOCUS.i, na.rm=T),
                       DEMI.PROFESS.i=mean(df$DEMI.PROFESS.i, na.rm=T),
                       DEMI.RESOURCES.i=mean(df$DEMI.RESOURCES.i, na.rm=T),
                       DEMI.EXPERIENCE.i=mean(df$DEMI.EXPERIENCE.i, na.rm=T),
                       DEMI.1=mean(df$DEMI.1, na.rm=T),
                       att_intl_bin_lag=mean(df$att_intl_bin_lag, na.rm=T),
                       ISO=medEff)
out2 <- merTools::predictInterval(M2, newdata=newdata, level=0.90,
                                  n.sims=1000, stat="mean", type = "probability", include.resid.var=FALSE, fix.intercept.variance=T)
out2 <- cbind(out2, newdata$rating.inv, newdata$rsf_pfi)
names(out2) <- c("fit","upr","lwr","rating.inv","rsf_pfi" )
out2$arena <- "Int'l"

# plot
out <- rbind(out1, out2)

ggplot(out, aes(x=rating.inv, y=fit, ymin=lwr, ymax= upr, shape=factor(rsf_pfi))) +
  geom_point(aes(shape = factor(rsf_pfi))) +
  geom_line(data=out, aes(linetype = factor(rsf_pfi))) +
  geom_ribbon(aes(ymin=lwr,ymax=upr), alpha=0.1) +
  facet_rep_grid( arena~. , switch = "both",  labeller = label_wrap_gen(width = 16, multi_line = TRUE))+
  labs(x = "Low   <-   Issue urgency   ->   High", y = "Predicted probability of attention") +
  scale_linetype_discrete(name  ="Press Freedom", breaks=c("0.9", "0.4"), labels=c("High","Low"))+
  scale_shape_discrete(name  ="Press Freedom", breaks=c("0.9", "0.4"), labels=c("High","Low"))+
  theme_bw()

ggsave(filename = "Fig 4.jpg", width = 4.5, height = 3.5, dpi = 600, units = "in")

rm(out)
rm(out1)
rm(out2)
rm(medEff)
rm(newdata)



# Table A.4: Pairwise Pearson correlation coefficients between dependent variables --------

t <- df[c("att_dom", "att_intl", "att_dom2", "att_dom_std", "att_intl_std", "att_dom2_std")]
rcorr(as.matrix(t))
rm(t)




# Figure A.1: Spread of survey responses by macro variables --------

dist <- df[!is.na(df$DEMI.EFFORT.i),]
dist <- dist[c("fhcategory",
               "PEItype",
               "unna_gdppc",
               "undp_hdi",
               "wdi_internet",
               "p_polity2",
               "unna_pop")]
dist$fhcategory <- ordered(dist$fhcategory ,c("Not Free","Partly Free","Free"))
dist$count <- 1
distfh <- dist[!is.na(dist$fhcategory),]
distpei <- dist[!is.na(dist$PEItype),]

p1 <- ggplot(distfh) + 
  stat_summary(aes(x=fhcategory, y = count), 
               fun.y = function(x) length(x) / length(unique(x)), 
               geom = "bar") +
  labs(x="Freedom House classification of respondent's country", y="Number of respondents") 

p2 <- ggplot(distpei) + 
  stat_summary(aes(x=PEItype, y = count), 
               fun.y = function(x) length(x) / length(unique(x)), 
               geom = "bar") +
  labs(x="PEI type of respondent's country", y="Number of respondents") 

p3 <- ggplot(dist, aes(unna_gdppc)) +
  geom_histogram(bins = 20) +
  labs(x="GDP per capita (PPP) of respondent's country", y="Number of respondents") 

p4 <- ggplot(dist, aes(undp_hdi)) +
  geom_histogram(bins = 20) +
  labs(x="Human Development Index (HDI) of respondent's country", y="Number of respondents") 

p5 <- ggplot(dist, aes(wdi_internet)) +
  geom_histogram(bins = 20) +
  labs(x="Internet penetration rate of respondent's country", y="Number of respondents") +
  theme_bw()

p6 <- ggplot(dist, aes(p_polity2)) +
  geom_histogram(bins = 17) +
  labs(x="Combined Polity IV Score of respondent's country", y="Number of respondents") 

ggsave(plot = multiplot(p1, p2, p6, p3, p4, p5, cols=2), filename = "Fig A1.jpg", width = 9, height = 10, dpi = 600, units = "in")

rm(dist)
rm(distfh)
rm(distpei)
rm(p1)
rm(p2)
rm(p3)
rm(p4)
rm(p5)
rm(p6)







# Figure A.2: Involvement of survey respondents in their organization --------

fam <- df[!is.na(df$DEMI.EFFORT.i),]
fam <- fam[c("A37_worked", "A38_strategic", "A39_operations")]

p1 <- ggplot(data=fam, aes(A37_worked)) + 
  geom_histogram(aes(x=A37_worked,y=..count..),binwidth=1, color="black", fill="white") +
  labs(x="'How long have you worked with your organization'", y="Number of respondents") +
  scale_x_continuous(breaks = c(1,2,3,4,5,6), labels = c("Less than a year",
                                                         "1-2 years",
                                                         "3-5 years",
                                                         "6-10 years",
                                                         "11-20 years",
                                                         "More than 20 years"))
summary(fam$A38_strategic)
skewness(fam$A38_strategic, na.rm = T)   #-2.556571
kurtosis(fam$A38_strategic, na.rm = T)   #10.42163
p2 <- ggplot(data=fam, aes(A38_strategic)) + 
  geom_histogram(aes(x=A38_strategic,y=..count..),binwidth=1, color="black", fill="white") +
  scale_x_continuous(labels = fmt_dcimals(0)) +
  annotate("text", x=2, y=100, label= "mean = 9.23     median = 10")  +
  annotate("text", x=2, y=80, label= "skewness = -2.557      kurtosis = 10.422") +
  labs(x="'How involved are you in your organization's strategic decisions' (0 = not al all involved; 10 = very involved)", y="Number of respondents")


summary(fam$A39_operations)
skewness(fam$A39_operations, na.rm = T)   #-2.098748
kurtosis(fam$A39_operations, na.rm = T)   #8.025871
p3 <- ggplot(data=fam, aes(A39_operations)) + 
  geom_histogram(aes(x=A39_operations,y=..count..),binwidth=1, color="black", fill="white") +
  scale_x_continuous(labels = fmt_dcimals(0)) +
  annotate("text", x=2, y=100, label= "mean = 9.06      median = 10")  +
  annotate("text", x=2, y=80, label= "skewness = -2.099      kurtosis = 8.026")  +
  labs(x="'How involved are you in your organization's operational activities' (0 = not al all involved; 10 = very involved)", y="Number of respondents")

ggsave(plot = multiplot(p1, p2, p3, cols=1), filename = "Fig A2.jpg", width = 9, height = 10, dpi = 600, units = "in")

rm(p1)
rm(p2)
rm(p3)
rm(fam)



# Figure A.4: Distribution of organizational independent variables --------

g1 <- ggplot(df, aes(DEMI.EFFORT.i)) +
  geom_histogram(binwidth=0.1) + 
  geom_histogram(aes(x = DEMI.EFFORT), binwidth=0.1, fill = "white", colour = "gray", alpha = 0.2) +
  annotate("text", x = 0.1, y = 95, label = "alpha = .75", hjust = 0) +
  annotate("text", x = 0.1, y = 85, label = "r = .98, p<.001", hjust = 0) +
  xlab("Media effort (0-1)") +
  ylab("No. of groups") +
  coord_cartesian(ylim = c(0, 100), xlim = c(0, 1)) +
  theme_minimal()

g2 <- ggplot(df, aes(DEMI.FOCUS.i)) +
  geom_histogram(binwidth=0.1) + 
  geom_histogram(aes(x = DEMI.FOCUS), binwidth=0.1, fill = "white", colour = "gray", alpha = 0.2) +
  annotate("text", x = 0.1, y = 95, label = "alpha = .83", hjust = 0) +
  annotate("text", x = 0.1, y = 85, label = "r = .66, p<.001", hjust = 0) +
  xlab("Specialization (0-1)") +
  ylab("No. of groups") +
  coord_cartesian(ylim = c(0, 100), xlim = c(0, 1)) +
  theme_minimal()

g3 <- ggplot(df, aes(DEMI.PROFESS.i)) +
  geom_histogram(binwidth=0.1) + 
  geom_histogram(aes(x = DEMI.PROFESS), binwidth=0.1, fill = "white", colour = "gray", alpha = 0.2) +
  annotate("text", x = 0.1, y = 95, label = "alpha = .67", hjust = 0) +
  annotate("text", x = 0.1, y = 85, label = "r = .72, p<.001", hjust = 0) +
  xlab("Professionalization (0-1)") +
  ylab("No. of groups") +
  coord_cartesian(ylim = c(0, 100), xlim = c(0, 1)) +
  theme_minimal()

ggsave(plot = multiplot(g1, g2, g3, cols=2), filename = "Fig A4.jpg", width = 6, height = 6, dpi = 600, units = "in")

rm(g1)
rm(g2)
rm(g3)




# Table A.7: Descriptive statistics for variables used in models M1-M4 --------

desc <- df[!is.na(df$DEMI.EFFORT.i),]
desc$att_dom_bin <- as.numeric(as.character(desc$att_dom_bin))
desc$att_intl_bin  <- as.numeric(as.character(desc$att_intl_bin))

desc <- desc[c("att_dom", "att_dom_std", "att_dom_bin",
             "att_intl", "att_intl_std", "att_intl_bin",
             "rsf_pfi", "rating.inv",
             "DEMI.EFFORT.i", "DEMI.FOCUS.i", "DEMI.PROFESS.i",
             "DEMI.RESOURCES.i", "DEMI.EXPERIENCE.i",
             "DEMI.1",
             "att_dom_bin_lag", "att_intl_bin_lag", "att_dom_lag", "att_intl_lag",
             "att_elect_dom", "att_elect_intl")]

stargazer(desc, 
          covariate.labels = c("Domestic news attention, raw count",
                               "Domestic news attention, standardized",
                               "Domestic news attention, dichotomized",
                               "Int'l news attention, raw count",
                               "Int'l news attention, standardized",
                               "Int'l news attention, dichotomized",
                               "Press freedom",
                               "Issue urgency (extent of electoral malpractice)",
                               "Media effort",
                               "Specialization",
                               "Professionalization",
                               "Resources",
                               "Longevity",
                               "Network membership",
                               "Lagged domestic news attention, dichotomized",
                               "Lagged int'l news attention, dichotomized",
                               "Lagged domestic news attention, raw count",
                               "Lagged int'l news attention, raw count",
                               "Domestic attention to elections (offset in M3)",
                               "Int'l attention to elections (offset in M4)"
          ),
          type = "text", out="Table A7.htm")

rm(desc)

# Table A.8:  Comparing alternative logistic regressions explaining the threshold of domestic news attention to groups  --------


htmlreg((list(M1null, M1a, M1b, M1c, M1d, M1e, M1)),
        single.row=TRUE,
        bold = 0.05,
        stars = c(0.001, 0.01, 0.05),
        leading.zero = F,
        inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        file="Table A8.doc")

rcompanion::nagelkerke(M1a, M1null)
rcompanion::nagelkerke(M1b, M1null)
rcompanion::nagelkerke(M1c, M1null)
rcompanion::nagelkerke(M1d, M1null)
rcompanion::nagelkerke(M1e, M1null)
rcompanion::nagelkerke(M1, M1null)


# Table A.9:  Comparing alternative logistic regressions explaining the threshold of international news attention to groups --------


htmlreg((list(M2null, M2a, M2b, M2c, M2d, M2e, M2)),
        single.row=TRUE,
        bold = 0.05,
        stars = c(0.001, 0.01, 0.05),
        leading.zero = F,
        inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
        head.tag = TRUE, body.tag = TRUE,
        file="Table A9.doc")

rcompanion::nagelkerke(M2a, M2null)
rcompanion::nagelkerke(M2b, M2null)
rcompanion::nagelkerke(M2c, M2null)
rcompanion::nagelkerke(M2d, M2null)
rcompanion::nagelkerke(M2e, M2null)
rcompanion::nagelkerke(M2, M2null)



# Table A.10: Comparing alternative count models explaining domestic news attention to groups --------


tab_model((list(M3p,M3hp,M3hnb, M3nb)),
          transform = NULL,
          order.terms = c(1,3,2,11,4,5,6,12,13,14,7,8,9,10),
          show.se = T,
          string.se = "SE",
          show.ci = F,
          show.aic = T,
          show.re.var = T,
          show.dev = T,
          show.obs = T,
          p.style = c("asterisk"),
          p.threshold = c(0.05, 0.01,0.001),
          file = "Table A10.html",
          use.viewer = F)

## Coefficients
finalcoefs3 <- data.frame(estimate = c(coef(summary(M3p))$cond[, "Estimate"],
                                      coef(summary(M3nb))$cond[, "Estimate"],
                                      coef(summary(M3hp))$cond[, "Estimate"],
                                      coef(summary(M3hnb))$cond[, "Estimate"]
),
se = c(coef(summary(M3p))$cond[, "Std. Error"],
       coef(summary(M3nb))$cond[, "Std. Error"],
       coef(summary(M3hp))$cond[, "Std. Error"],
       coef(summary(M3hnb))$cond[, "Std. Error"]
),
pvalue = c(coef(summary(M3p))$cond[, "Pr(>|z|)"],
           coef(summary(M3nb))$cond[, "Pr(>|z|)"],
           coef(summary(M3hp))$cond[, "Pr(>|z|)"],
           coef(summary(M3hnb))$cond[, "Pr(>|z|)"]
),
model    = rep(c("Pois", "NB", "Hurdle-Pois", "Hurdle-NB"), each = 14),
term     = rep(names(coef(summary(M3p))$cond[, "Estimate"]), 4))



zerocoefs3 <- data.frame(estimate = c(coef(summary(M3hp))$zi[, "Estimate"],
                                     coef(summary(M3hnb))$zi[, "Estimate"]
),
se = c(coef(summary(M3hp))$zi[, "Std. Error"],
       coef(summary(M3hnb))$zi[, "Std. Error"]
),
pvalue = c(coef(summary(M3hp))$zi[, "Pr(>|z|)"],
           coef(summary(M3hnb))$zi[, "Pr(>|z|)"]
),
model    = rep(c("Hurdle-Pois", "Hurdle-NB"), each = 7),
term     = rep(names(coef(summary(M3zinb))$zi[, "Estimate"]), 2))


newdata <- dfc3[complete.cases(dfc3[c("att_dom",
                                              "DEMI.RESOURCES.i",
                                              "DEMI.EFFORT.i",
                                              "DEMI.FOCUS.i",
                                              "DEMI.PROFESS.i",
                                              "DEMI.1",
                                              "rsf_pfi",
                                              "rating.inv",
                                              "att_dom_lag",
                                              'att_elect_dom')]),]

pred <- predict(M3p, newdata, se.fit=T)
pred$M3p<- pred$fit
newdata <- cbind(newdata, pred[3])

pred <- predict(M3nb, newdata, se.fit=T)
pred$M3nb <- pred$fit
newdata <- cbind(newdata, pred[3])

pred <- predict(M3hp, newdata, se.fit=T)
pred$M3hp <- pred$fit
newdata <- cbind(newdata, pred[3])

pred <- predict(M3hnb, newdata, se.fit=T)
pred$M3hnb <- pred$fit
newdata <- cbind(newdata, pred[3])

#  root mean squared error
rmse(newdata$att_dom-newdata$M3p) 
rmse(newdata$att_dom-newdata$M3hp) 
rmse(newdata$att_dom-newdata$M3hnb)  
rmse(newdata$att_dom-newdata$M3nb)

rm(dfc3)
rm(newdata)
rm(finalcoefs3)
rm(zerocoefs3)

# Table A.11: Comparing alternative count models explaining international news attention to groups --------

tab_model((list(M4p,M4hp,M4hnb, M4nb)),
          transform = NULL,
          order.terms = c(1,3,2,11,4,5,6,12,13,14,7,8,9,10),
          show.se = T,
          string.se = "SE",
          show.ci = F,
          show.aic = T,
          show.re.var = T,
          show.dev = T,
          show.obs = T,
          p.style = c("asterisk"),
          p.threshold = c(0.05, 0.01,0.001),
          file = "Table A11.html",
          use.viewer = F)

# coefficients
finalcoefs4 <- data.frame(estimate = c(coef(summary(M4p))$cond[, "Estimate"],
                                      coef(summary(M4nb))$cond[, "Estimate"],
                                      coef(summary(M4hp))$cond[, "Estimate"],
                                      coef(summary(M4hnb))$cond[, "Estimate"]
),
se = c(coef(summary(M4p))$cond[, "Std. Error"],
       coef(summary(M4nb))$cond[, "Std. Error"],
       coef(summary(M4hp))$cond[, "Std. Error"],
       coef(summary(M4hnb))$cond[, "Std. Error"]
),
pvalue = c(coef(summary(M4p))$cond[, "Pr(>|z|)"],
           coef(summary(M4nb))$cond[, "Pr(>|z|)"],
           coef(summary(M4hp))$cond[, "Pr(>|z|)"],
           coef(summary(M4hnb))$cond[, "Pr(>|z|)"]
),
model    = rep(c("Pois", "NB", "Hurdle-Pois", "Hurdle-NB"), each = 14),
term     = rep(names(coef(summary(M4p))$cond[, "Estimate"]), 4))



zerocoefs4 <- data.frame(estimate = c(coef(summary(M4hp))$zi[, "Estimate"],
                                     coef(summary(M4hnb))$zi[, "Estimate"]
),
se = c(coef(summary(M4hp))$zi[, "Std. Error"],
       coef(summary(M4hnb))$zi[, "Std. Error"]
),
pvalue = c(coef(summary(M4hp))$zi[, "Pr(>|z|)"],
           coef(summary(M4hnb))$zi[, "Pr(>|z|)"]
),
model    = rep(c("Hurdle-Pois", "Hurdle-NB"), each = 7),
term     = rep(names(coef(summary(M4hnb))$zi[, "Estimate"]), 2))


#  root mean squared error
newdata <- dfc4[complete.cases(dfc4[c("att_intl",
                                      "DEMI.RESOURCES.i",
                                      "DEMI.EFFORT.i",
                                      "DEMI.FOCUS.i",
                                      "DEMI.PROFESS.i",
                                      "DEMI.1",
                                      "rsf_pfi",
                                      "rating.inv",
                                      "att_intl_lag",
                                      "att_elect_intl")]),]

pred <- predict(M4p, newdata, se.fit=T)
pred$M4p<- pred$fit
newdata <- cbind(newdata, pred[3])
pred <- predict(M4nb, newdata, se.fit=T)
pred$M4nb <- pred$fit
newdata <- cbind(newdata, pred[3])
pred <- predict(M4hp, newdata, se.fit=T)
pred$M4hp <- pred$fit
newdata <- cbind(newdata, pred[3])
pred <- predict(M4hnb, newdata, se.fit=T)
pred$M4hnb <- pred$fit
newdata <- cbind(newdata, pred[3])

rmse(newdata$att_intl-newdata$M4p)    
rmse(newdata$att_intl-newdata$M4hp)  
rmse(newdata$att_intl-newdata$M4hnb) 
rmse(newdata$att_intl-newdata$M4nb)

rm(finalcoefs4)
rm(zerocoefs4)
rm(newdata)
rm(dfc4)
rm(pred)




# Table A.12: Models without lagged attention --------


M1z <- glmer(att_dom_bin ~ scale(rsf_pfi)*scale(rating.inv) +
              scale(rsf_pfi)*scale(DEMI.EFFORT.i)  + scale(rsf_pfi)*scale(DEMI.FOCUS.i) + scale(rsf_pfi)*scale(DEMI.PROFESS.i) + 
              scale(DEMI.RESOURCES.i) +scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) +
              (1 | ISO), data=df, family=binomial(link = "logit"),
            control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))


M2z <- glmer(att_intl_bin ~ scale(rsf_pfi)*scale(rating.inv) +
              scale(rsf_pfi)*scale(DEMI.EFFORT.i)  + scale(rsf_pfi)*scale(DEMI.FOCUS.i) + scale(rsf_pfi)*scale(DEMI.PROFESS.i) + 
              scale(DEMI.RESOURCES.i) +scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) +
              (1 | ISO), data=df, family=binomial(link = "logit"),
            control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))


dfc3 <- df[complete.cases(df[c("att_dom",
                               "DEMI.RESOURCES.i",
                               "DEMI.EFFORT.i",
                               "DEMI.FOCUS.i",
                               "DEMI.PROFESS.i",
                               "DEMI.EXPERIENCE.i",
                               "DEMI.1",
                               "rsf_pfi",
                               "rating.inv",
                               "att_dom_lag",
                               "p_polity2",
                               "att_elect_dom",
                               "unna_gdppc",
                               "p_polity2",
                               "ISO")]),]

form <- att_dom ~ scale(rsf_pfi)*scale(rating.inv) + 
  scale(rsf_pfi)*scale(DEMI.EFFORT.i)  + scale(rsf_pfi)*scale(DEMI.FOCUS.i) + scale(rsf_pfi)*scale(DEMI.PROFESS.i) + 
  scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) + 
  offset(log(att_elect_dom)) +
  (1|ISO)

M3nbz = glmmTMB(form, data=dfc3, family=nbinom2, control = glmmTMBControl(optCtrl=list(iter.max=1e9,eval.max=1e9)))


dfc4 <- df[complete.cases(df[c("att_intl",
                               "DEMI.RESOURCES.i",
                               "DEMI.EFFORT.i",
                               "DEMI.FOCUS.i",
                               "DEMI.PROFESS.i",
                               "DEMI.EXPERIENCE.i",
                               "DEMI.1",
                               "rsf_pfi",
                               "rating.inv",
                               "att_intl_lag",
                               "att_elect_intl",
                               "unna_gdppc",
                               "p_polity2",
                               "ISO")]),]

form <- att_intl ~ scale(rsf_pfi)*scale(rating.inv) + 
  scale(rsf_pfi)*scale(DEMI.EFFORT.i)  + scale(rsf_pfi)*scale(DEMI.FOCUS.i) + scale(rsf_pfi)*scale(DEMI.PROFESS.i) + 
  scale(DEMI.RESOURCES.i) + scale(DEMI.EXPERIENCE.i) + scale(DEMI.1) +
  offset(log(att_elect_intl)) +
  (1|ISO)

M4nbz = glmmTMB(form, data=dfc4, family=nbinom2, control = glmmTMBControl(optCtrl=list(iter.max=1e9,eval.max=1e9)))



tab_model((list(M1z, M2z, M3nbz, M4nbz)),
          transform = NULL,
          show.se = T,
          string.se = "SE",
          show.ci = F,
          show.aic = T,
          show.re.var = T,
          show.dev = T,
          show.obs = T,
          p.style = c("asterisk"),
          p.threshold = c(0.05, 0.01,0.001),
          file = "Table A12.html",
          use.viewer = F)
